Here we show the code to reproduce the analyses of: Risso and Pagnotta (2020). Per-sample standardization and asymmetric winsorization lead to accurate classification of RNA-seq expression profiles.
This file belongs to the repository: https://github.com/drisso/awst_analysis.
The code is released with license GPL v3.0.
Install and load awst
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("drisso/awst")
rm(list = ls())
library(readr)
library(awst)
library(EDASeq)
library(Rtsne)
library(umap)
library(dendextend)
setwd("~/Dropbox/AWST/CITE20201018/")
jobName <- "CBMC20201018"
Data import and cleaning
The data are available on GEO with the following accession number.
In the following chunks of code, we arrange a data-frame where each of the cell are annotated, and the annotations are coded as colors.
We mostly follow the original article for data preprocessing.
Note that here we download and read the data directly from the remote gzipped files available in GEO. In some cases, with unstable internet connections, the following chunks may not work. In such cases, we suggest to “manually” downaload the data using the links provided above and read the data in R from a local copy of the file.
ADT and annotation
if(!file.exists("annotation.RData")) {
positive.col <- "red"
negative.col <- "gray"
ttable <- read_csv("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE100866&format=file&file=GSE100866_PBMC_vs_flow_10X-ADT_umi.csv.gz")
ttable <- as.data.frame(ttable)
rownames(ttable) <- ttable[,1]
ttable <- ttable[,-1]
ttable <- 1 + t(ttable)
tmp <- exp(log(apply(ttable, 1, prod))/ncol(ttable))
ttable <- log(ttable/tmp)
ttable <- 6 * (pnorm(ttable) - 0.5)
annotation.df <- as.data.frame(ttable)
annotation.df$cell <- rownames(annotation.df)
nBins <- 9
bbreaks <- seq(-3, 3, length.out = nBins + 1)
levels_cols <- c("green4", "green3", "green2", "green", "gray75", "red", "red2", "red3", "red4")
j <- 0
while(j < ncol(ttable)) {
j <- j + 1
annotation.df$tmp <- cut(ttable[, j], breaks = bbreaks, include.lowest = TRUE)
levels(annotation.df$tmp) <- levels_cols
colnames(annotation.df)[ncol(annotation.df)] <- paste0(colnames(ttable)[j], ".col")
}
save(annotation.df, file = "annotation.RData")
} else load("annotation.RData")
Single-cell RNA-seq
Before applying normalization and AWST, we filtered out cells that did not met the following criteria:
if(!file.exists("Level3.RData")) {
ddata <- read_csv("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE100866&format=file&file=GSE100866_PBMC_vs_flow_10X-RNA_umi.csv.gz")
#restrict the features to the HUMAN ones
ddata <- ddata[grep("^HUMAN_", ddata$X1),]
ddata <- as.data.frame(ddata)
rownames(ddata) <- gsub("HUMAN_", "", ddata$X1)
ddata <- as.matrix(ddata[,-1])
save(ddata, file = "Level3.RData")
} else if(!file.exists("normCounts.RData")) load("Level3.RData")
Apply AWST
if(!file.exists("normCounts.RData")) {
dim(ddata <- t(ddata)) # 7985 17014
###
no_of_detected_gene_per_sample <- rowSums(ddata > 0)
fivenum(no_of_detected_gene_per_sample)# 10 638 739 873 4833
# restrict the collection of cells to those cells having at least 500 observed features
sum(no_of_detected_gene_per_sample > 500) # 7613
dim(ddata <- ddata[no_of_detected_gene_per_sample > 500,])#[1] 7613 17014
# apply the full quantile normalization
normCounts <- EDASeq::betweenLaneNormalization(t(ddata), which = "full", round = FALSE)
save(normCounts, file = "normCounts.RData")
} else if(!file.exists("expression.RData")) load("normCounts.RData")
# apply the AWS-transformation
if(!file.exists("expression.RData")) {
exprData <- awst(normCounts, poscount = TRUE, full_quantile = TRUE)
dim(exprData <- gene_filter(exprData)) #[1] 7613 230
save(exprData, file = "expression.RData")
} else load("expression.RData")
Clustering and dimensionality reduction
Note that due to changes to the pseudo-random number generator in R 3.6, the behavior of set.seed() has changed. Hence, the t-SNE and UMAP plots below are not exact copies of the ones in the paper, obtained with an older version of R. However, the main features of the datasets are preserved.
if(!file.exists("expression_dist_hclust.RData")) {
nrow_exprData <- nrow(exprData)
ncol_exprData <- ncol(exprData)
ddist <- dist(exprData)
# save(ddist, nrow_exprData, ncol_exprData, file = "expression_dist.RData")
hhc <- hclust(ddist, method = "ward.D2")
save(hhc, nrow_exprData, ncol_exprData, file = "expression_dist_hclust.RData")
} else load("expression_dist_hclust.RData")
if(!file.exists("expression_prcomp.RData")) {
pprcomp <- prcomp(exprData)
pprcomp$x <- pprcomp$x[, 1:10]
pprcomp$rotation <- pprcomp$rotation[, 1:10]
save(pprcomp, file = "expression_prcomp.RData")
} else load("expression_prcomp.RData")
if(!file.exists("expression_Rtsne_2d.RData")) {
set.seed(2019) # needed to get the figure in the paper
ans_Rtsne <- Rtsne(exprData, pca = FALSE) # Run TSNE
save(ans_Rtsne, file = "expression_Rtsne_2d.RData")
} else load("expression_Rtsne_2d.RData")
if(!file.exists("expression_umap_2d.RData")) {
set.seed(2019) # needed to get the figure in the paper
ans_umap <- umap(exprData)
save(ans_umap, file = "expression_umap_2d.RData")
} else load("expression_umap_2d.RData")
load("annotation.RData")
load("expression_dist_hclust.RData")
annotation.df <- annotation.df[hhc$labels,]
###############
save_plots <- FALSE
png_width_large <- 2100
png_height_large <- 500
png_width_small <- width_png <- 700
png_height_small <- 700
png_res <- 1/300
###################
color.bar2 <- function(x_pos, y_pos, lut, min, max=-min, nticks=11, ticks=seq(min, max, len=nticks), title='', values = NULL) {
scale = (max-min)/length(lut)*0.3
for (i in 1:length(lut)) {
y_low <- (i-1)*scale + min + y_pos
y_high <- y_low + scale
rect(x_pos,y_low,x_pos+.05,y_high, col=lut[i], border=NA)
text(x_pos+.05, (y_low + y_high)/2, values[i], adj = -0.1)
}
}
color.bar3 <- function(x_pos, y_pos, lut, min, max=-min, nticks=11, ticks=seq(min, max, len=nticks), title='', values = NULL) {
scale = (max-min)/length(lut)*0.3
for (i in 1:length(lut)) {
y_low <- (i-1)*scale + y_pos
y_high <- y_low + scale
rect(x_pos,y_low,x_pos+90,y_high, col=lut[i], border=NA)
text(x_pos+.05, (y_low + y_high)/2, values[i], adj = -2)
}
}
vvalues <- c("-3.0", "-2.0", "-1.3", "-0.6", " 0.0", " 0.6", " 1.3", " 2.0", " 3.0")
ffill2 <- names(table(annotation.df$CD3.col))
Main Clustering
clustering.prefix <- "CBMC"; short.prefix <- "CBMC"
clustering.df <- data.frame(cell = annotation.df$cell)
rownames(clustering.df) <- clustering.df$cell
#clustering.df <- clustering.df[hhc$labels,]
############
#load(paste0(jobName, "_dist_hclust.RData"))
mmain <- paste0("(a) - CBMC study (", nrow_exprData, " cells/", ncol_exprData, " genes)")
if(save_plots) {
# mmain <- ""
png(paste0(jobName, "_expression_dist_hclust.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
#mmain = paste0(jobName, " (", nrow_exprData, " cells/", ncol_exprData, " genes)")
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
abline(h = hh, col = "red")
tmp <- tmp_ <- as.factor(cutree(hhc, k = wwhere))
wwhere <- length(unique(levels(tmp)))
clusteringWhere <- paste0(clustering.prefix, wwhere)
clusteringWhere.col <- paste0(clusteringWhere, ".col")
assign(clusteringWhere.col, tmp)
levels(tmp) <- paste0(short.prefix, wwhere, 1:wwhere)
assign(clusteringWhere, tmp)
levels(tmp) <- palette()[1:wwhere]
assign(clusteringWhere.col, tmp)
tt <- table(get(clusteringWhere), get(clusteringWhere.col))
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
names(colorCode) <- rownames(tt)
assign(paste0(clusteringWhere, ".colorCode"), colorCode)
clust.colorCode <- colorCode
clustering.df$tmp <- get(clusteringWhere)
clustering.df$tmp.explanatory <- clustering.df$tmp
clustering.df$tmp.col <- get(clusteringWhere.col)
ncol_ <- ncol(clustering.df)
colnames(clustering.df)[(ncol_-2):ncol_] <- c(clusteringWhere, paste0(clusteringWhere, ".explanatory"), clusteringWhere.col)
levels(clustering.df[, paste0(clusteringWhere, ".explanatory")]) <- c("CBMC1 - T Cell", "CBMC2 - B Cell", "CBMC3 - unclear", "CBMC4 - Natutal Killer", "CBMC5 - Monocyte", "CBMC6 - myeloid DC")
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- clustering.df[, ncol(clustering.df)]
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
save(clustering.df, clust.colorCode, file = "clustering.RData")
tt <- table(clustering.df[, ncol(clustering.df)-1])
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- AWST_legend <- paste(names(tt), " (", tt, "; ", pct, ")", sep = "")
tt <- table(clustering.df[, ncol(clustering.df)-1], clustering.df[, ncol(clustering.df)])
ffill <- AWST_fill <- colnames(tt)[apply(tt, 1, which.max)]
legend(1, .85, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
#text(1, 1, "(a)")
####
#fName <- paste0(jobName, "_metadata.tsv")
#write.table(annotation.col, file = fName, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
PCA (AWST) | …
x_legend <- 0.08; y_legend <- -1.75
x_text <- -2.3; y_text <- 2.
load("expression_prcomp.RData")
pprcomp$x <- scale(pprcomp$x)
#
fName <- paste0(jobName, "_expression_prcomp.tsv")
write.table(pprcomp$x[, 1:3], file = fName, sep = "\t", row.names = FALSE, col.names = FALSE)
#
mmain = paste0("Principal components analysis (", nrow_exprData, " cells/", ncol_exprData, " features)")
cat(sprintf("\n\n### %s\n\n", mmain))
Principal components analysis (7613 cells/230 features)
if(save_plots) png(paste0(jobName, "_expression_prcomp_CITE6.png"), width= png_width_small, height= png_height_small, res = png_res)
plot(pprcomp$x, col = clustering.df$CBMC6.col, main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
legend(x_legend, y_legend, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
text(x_text, y_text, "(b)")
###
mmain <- "Principal component analysis - CD3"# "prcomp (AWST)| flow cytometry/CD3"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD3
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD3.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD3.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
ffill2 <- names(table(annotation.df$CD3.col))
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(a)")
#
mmain <- "Principal component analysis - CD19"# "prcomp (AWST)| flow cytometry/CD19"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD19
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD19.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD19.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(b)")
#
mmain <- "Principal component analysis - CD11c"# "prcomp (AWST)| flow cytometry/CD11c"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD11c
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD11c.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD11c.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(c)")
#
mmain <- "Principal component analysis - CD14"# "prcomp (AWST)| flow cytometry/CD14"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD14
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD14.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD14.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(d)")
######### CD4
mmain <- "Principal component analysis - CD4"# "prcomp (AWST)| flow cytometry/CD4"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD4
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD4.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD4.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(e)")
######### CD8
mmain <- "Principal component analysis - CD8"# "prcomp (AWST)| flow cytometry/CD8"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD8
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD8.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD8.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(f)")
######### CD2
mmain <- "Principal component analysis - CD2"# "prcomp (AWST)| flow cytometry/CD8"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD2
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD2.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD2.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(g)")
######### CD57
mmain <- "Principal component analysis - CD57"# "prcomp (AWST)| flow cytometry/CD8"
cat(sprintf("\n\n### %s\n\n", mmain))
Principal component analysis - CD57
if(save_plots) png(paste0(jobName, "_expression_prcomp_CD57.png"), width = width_png, height = width_png)
plot(pprcomp$x, col = as.character(annotation.df$CD57.col), main = mmain,
xlab = "first principal component", ylab = "second principal component", pch = 19)
color.bar2(-2.3, 0.3, ffill2, -3, values = vvalues)
text(x_text, y_text, "(h)")
Rtsne (AWST) | clustering
load("expression_Rtsne_2d.RData")
ans_Rtsne$Y <- scale(ans_Rtsne$Y)
if(save_plots) png(paste0(jobName, "_expression_Rtsne.png"), width= png_width_small, height= png_height_small, res = png_res)
mmain = "T-distributed Stochastic Neighbor Embedding (tsne)"
plot(-ans_Rtsne$Y[, 1], ans_Rtsne$Y[, 2], col = clustering.df$CBMC6.col, main = mmain, xlab = "", ylab = "")
legend(-2.5, -1.1, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
text(-2.5, 1.8, "(c)")
umap (AWST) | clustering
load("expression_umap_2d.RData")
ans_umap$layout <- scale(ans_umap$layout)
if(save_plots) png(paste0(jobName, "_expression_umap.png"), width= png_width_small, height= png_height_small, res = png_res)
mmain = "Uniform Manifold Approximation and Projection (umap)"
plot(-ans_umap$layout[, 1], -ans_umap$layout[, 2], col = clustering.df$CBMC6.col, main = mmain, xlab = "", ylab = "")
legend(1, -1.5, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
text(-1.5, 1, "(d)")
Compute scran HVG
if(!file.exists("scran_hvg.RData")) {
load("Level3.RData")
library(scran)
no_of_detected_gene_per_sample <- colSums(ddata > 0)
ddata <- ddata[,no_of_detected_gene_per_sample > 500]
sfs <- scran::calculateSumFactors(ddata)
lognorm <- scater::normalizeCounts(ddata, size_factors = sfs)
stats <- modelGeneVar(lognorm)
scran_hvg <- getTopHVGs(stats, fdr.threshold=0.05)
save(scran_hvg, file = "scran_hvg.RData")
} else load("scran_hvg.RData")
AWST with scran
if(!file.exists("AWST_expression_dist_hclust_scran.RData")) {
load("normCounts.RData")
awstData <- awst(normCounts, poscount = TRUE, full_quantile = TRUE)
awst_scran <- awstData[, scran_hvg]
nrow_awstData <- nrow(awst_scran)
ncol_awstData <- ncol(awst_scran)
ddist <- dist(awst_scran)
hhc <- hclust(ddist, method = "ward.D2")
save(hhc, nrow_awstData, ncol_awstData, scran_hvg, file = "AWST_expression_dist_hclust_scran.RData")
} else load("AWST_expression_dist_hclust_scran.RData")
mmain <- paste0("(b) - CBMC study (", nrow_awstData, " cells/", ncol_awstData, " genes), AWST & scran")
if(save_plots) {
# mmain <- ""
png(paste0(jobName, "_expression_dist_hclust_AWST_scran.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'AWST/scran' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, clust.colorCode, file = "clustering_awst_scran.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
llegend <- llegend[AWST_fill]
ffill <- ffill[AWST_fill]
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "AWST & scran", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Analysis with VST normalization
if(!file.exists("VST_expression.RData")) {
library(DESeq2)
load("Level3.RData")
no_of_detected_gene_per_sample <- colSums(ddata > 0)
ddata <- ddata[,no_of_detected_gene_per_sample > 500]
# ddata[ddata == 0] <- 1
annotation.df <- annotation.df[colnames(ddata),]
dds <- DESeqDataSetFromMatrix(countData = as.matrix(ddata),
colData = annotation.df,
design = ~1)
dds <- estimateSizeFactors(dds, type="poscounts")
# system.time({vvst <- vst(dds, blind=FALSE)})
# Error in vst(dds, blind = FALSE) :
# less than 'nsub' rows with mean normalized count > 5,
# it is recommended to use varianceStabilizingTransformation directly
system.time({vvst <- varianceStabilizingTransformation(dds, blind=FALSE)})
#-- note: fitType='parametric', but the dispersion trend was not well captured by the
# function: y = a/x + b, and a local regression fit was automatically substituted.
# specify fitType='local' or 'mean' to avoid this message next time.
# user system elapsed
#1006.131 20.738 1009.941
save(vvst, file = "VST_expression.RData")
} else load("VST_expression.RData")
if(!file.exists("VST_expression_dist_hclust.RData")) {
vars <- matrixStats::rowVars(assay(vvst))
names(vars) <- rownames(vvst)
vars <- sort(vars, decreasing = TRUE)
vst_hvg <- assay(vvst)[names(vars)[seq_len(ncol(exprData))],]
nrow_vst_hvg <- ncol(vst_hvg)
ncol_vst_hvg <- nrow(vst_hvg)
dist_vst_hvg <- dist(t(vst_hvg))
hhc_hvg <- hclust(dist_vst_hvg, method = "ward.D2")
vst_gf <- assay(vvst)[colnames(exprData),]
nrow_vst_gf <- ncol(vst_gf)
ncol_vst_gf <- nrow(vst_gf)
dist_vst_gf <- dist(t(vst_gf))
hhc_gf <- hclust(dist_vst_gf, method = "ward.D2")
vst_scran <- assay(vvst)[scran_hvg,]
nrow_vst_scran <- ncol(vst_scran)
ncol_vst_scran <- nrow(vst_scran)
dist_vst_scran <- dist(t(vst_scran))
hhc_scran <- hclust(dist_vst_scran, method = "ward.D2")
save(nrow_vst_hvg, ncol_vst_hvg, hhc_hvg,
nrow_vst_gf, ncol_vst_gf, hhc_gf,
nrow_vst_scran, ncol_vst_scran, hhc_scran,
file = "VST_expression_dist_hclust.RData")
} else load("VST_expression_dist_hclust.RData")
VST & Highly Variable Genes
hhc <- hhc_hvg
mmain <- paste0("(a) - CBMC study (", nrow_vst_gf, " cells/", ncol_vst_gf, " genes), VST & HVG")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_VST_HVG.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[4] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'VST/HVG' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_VST_HVG.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "VST & HVG", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
VST & genes selected by AWST
hhc <- hhc_gf
mmain <- paste0("(b) - CBMC study (", nrow_vst_gf, " cells/", ncol_vst_gf, " genes), VST & genes from AWST")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_VST_GF.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[4] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'VST/HVG from AWST' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_VST_GF.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "VST/HVG from AWST", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
VST & scran
hhc <- hhc_scran
mmain <- paste0("(c) - CBMC study (", nrow_vst_gf, " cells/", ncol_vst_gf, " genes), VST & scran")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_VST_scran.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[6] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'VST/scran' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_VST_scran.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "VST/scran", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Analysis with Townes normalization
if(!file.exists("SCRY_expression.RData")) {
library(scry)
load("Level3.RData")
no_of_detected_gene_per_sample <- colSums(ddata > 0)
ddata <- ddata[,no_of_detected_gene_per_sample > 500]
# ddata[ddata == 0] <- 1
system.time(scryData_pearson <- nullResiduals(as.matrix(ddata), type = "pearson"))
# user system elapsed
# 5.909 3.773 4.441
ans_devianceFeatureSelection <- devianceFeatureSelection(as.matrix(ddata))
names(ans_devianceFeatureSelection) <- rownames(ddata)
ans_devianceFeatureSelection <- sort(ans_devianceFeatureSelection, decreasing = TRUE)
save(scryData_pearson, ans_devianceFeatureSelection, file = "SCRY_expression.RData")
} else load("SCRY_expression.RData")
if(!file.exists("SCRY_expression_dist_hclust.RData")) {
# select high variable genes
vars <- matrixStats::rowVars(scryData_pearson)
names(vars) <- rownames(scryData_pearson)
vars <- sort(vars, decreasing = TRUE)
scry_hvg <- scryData_pearson[names(vars)[seq_len(ncol(exprData))],]
nrow_hvg <- ncol(scry_hvg)
ncol_hvg <- nrow(scry_hvg)
dist_hvg <- dist(t(scry_hvg))
hhc_hvg <- hclust(dist_gf, method = "ward.D2")
scry_gf <- scryData_pearson[colnames(exprData),]
nrow_gf <- ncol(scry_gf)
ncol_gf <- nrow(scry_gf)
dist_gf <- dist(t(scry_gf))
hhc_gf <- hclust(dist_gf, method = "ward.D2")
scry_scran <- scryData_pearson[scran_hvg,]
nrow_scran <- ncol(scry_scran)
ncol_scran <- nrow(scry_scran)
dist_scran <- dist(t(scry_scran))
hhc_scran <- hclust(dist_scran, method = "ward.D2")
hvg_scry <- names(ans_devianceFeatureSelection[seq_len(ncol(exprData))])
scry_scry <- scryData_pearson[hvg_scry,]
nrow_scry <- ncol(scry_scry)
ncol_scry <- nrow(scry_scry)
dist_scry <- dist(t(scry_scry))
hhc_scry <- hclust(dist_scry, method = "ward.D2")
save(nrow_hvg, ncol_hvg, hhc_hvg,
nrow_gf, ncol_gf, hhc_gf,
nrow_scran, ncol_scran, hhc_scran,
hhc_scry, nrow_scry, ncol_scry, hvg_scry,
file = "SCRY_expression_dist_hclust.RData")
} else load("SCRY_expression_dist_hclust.RData")
Townes & Highly Variable Genes
hhc <- hhc_hvg
mmain <- paste0("(a) - CBMC study (", nrow_gf, " cells/", ncol_gf, " genes), Townes & HVG")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_scry_HVG.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[c(2, 4:6)] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'Townes/HVG' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_scry_HVG.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
#llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
#ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "VST & HVG", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Townes & genes selected by AWST
hhc <- hhc_gf
mmain <- paste0("(b) - CBMC study (", nrow_gf, " cells/", ncol_gf, " genes), Townes & genes from AWST")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_scry_GF.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[5:6] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'Townes/HVG from AWST' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_scry_GF.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
llegend[6] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
ffill[6] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "Townes/HVG from AWST", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Townes & scran
hhc <- hhc_scran
mmain <- paste0("(c) - CBMC study (", nrow_scran, " cells/", ncol_scran, " genes), Townes & scran")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_scry_scran.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[c(2, 4:6)] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'Townes/scran' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
#save(clust.df, colorCode, file = "clustering_scry_scran.RData")
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
#llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
#ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "Townes/scran", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Townes & scry
hhc <- hhc_scry
mmain <- paste0("(d) - CBMC study (", nrow_scran, " cells/", ncol_scran, " genes), Townes & scry")
if(save_plots) {
png(paste0(jobName, "_expression_dist_hclust_scry_dev.png"),
width= png_width_large, height= png_height_large, res = png_res)
}
hhc$height <- hhc$height/max(hhc$height)
annotation.df <- annotation.df[hhc$labels,]
clust.df <- data.frame(cell = clustering.df$cell,
CBMC6.explanatory = clustering.df$CBMC6.explanatory,
CBMC6.col = clustering.df$CBMC6.col,
row.names = clustering.df$cell)
plot(hhc, hang = -1, labels = FALSE, xlab = "", sub = "", main = mmain)
###
wwhere <- 6
hh <- mean(c(hhc$height[length(hhc$height)-wwhere+2], hhc$height[length(hhc$height)-wwhere+1]))
segments(1, hh, nrow_exprData, hh, col = "red")
clust.df$tmp <- as.factor(cutree(hhc, k = wwhere))
tt <- table(clust.df$tmp, clust.df$CBMC6.col)
colorCode <- colnames(tt)[apply(tt, 1, which.max)]
colorCode[5:6] <- "gray"
levels(clust.df$tmp) <- colorCode
annotation.col <- annotation.df[, grep(".col", colnames(annotation.df))]
colnames(annotation.col) <- gsub(".col", "", colnames(annotation.col))
annotation.col <- annotation.col[, rev(c( "CD19", "CD3", "CD11c", "CD14","CD4","CD8", "CD2","CD57"))]
annotation.col$CBMC <- as.character(clust.df$CBMC6.col)
annotation.col$'Townes/scry' <- as.character(clust.df$tmp)
colored_bars(colors = annotation.col, dend = as.dendrogram(hhc), y_scale = 0.29, y_shift = 0.0045)
tt <- table(clust.df$tmp)
pct <- paste0(round(100*tt/sum(tt), 1), "%")
llegend <- paste("", tt, " cells; (", pct, ")", sep = "")
ffill <- names(tt)
names(llegend) <- names(ffill) <- ffill
tmp <- llegend
llegend <- llegend[AWST_fill]
#llegend[is.na(llegend)] <- setdiff(tmp, llegend)
tmp <- ffill
ffill <- ffill[AWST_fill]
#ffill[is.na(ffill)] <- setdiff(tmp, ffill)
legend(1, 0.97, legend=AWST_legend, fill = AWST_fill, y.intersp = 1, box.col = "white", border = "white", title = "CBMC", title.adj = 0)
legend(1500, 0.97, legend=llegend, fill = ffill, y.intersp = 1, box.col = "white", border = "white", title = "Townes/scry", title.adj = 0)
color.bar3(7230, 0.65, ffill2, -0.6, values = vvalues)
text(7275, 0.63, "Marker's level")
Feature selection steps
Here we compare the genes selected by each approach.
library(UpSetR)
# vars <- matrixStats::rowVars(ddata)
# names(vars) <- rownames(ddata)
# vars <- sort(vars, decreasing = TRUE)
load("expression.RData")
forUpset <- list("AWST +\n geneFilter" = colnames(exprData),
"scran HVG" = scran_hvg,
"scry HVG" = hvg_scry)
if(save_plots) {
png(paste0(jobName, "_upset.png"),
width= 1100, height= 700, res = png_res)
}
upset(fromList(forUpset), order.by = "freq",
main.bar.color = "#1F78B4", sets.bar.color = "#1F78B4")
library(eulerr)
if(save_plots) {
png(paste0(jobName, "_eulerr.png"),
width= png_height_large, height= png_height_large, res = png_res)
}
plot(euler(forUpset),quantities = TRUE)
#knitr::knit_exit()
Session info
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=it_IT.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=it_IT.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] eulerr_6.1.0 UpSetR_1.4.0
## [3] dendextend_1.13.4 umap_0.2.5.0
## [5] Rtsne_0.15 EDASeq_2.18.0
## [7] ShortRead_1.42.0 GenomicAlignments_1.20.1
## [9] SummarizedExperiment_1.14.1 DelayedArray_0.10.0
## [11] matrixStats_0.56.0 Rsamtools_2.0.3
## [13] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0
## [15] Biostrings_2.52.0 XVector_0.24.0
## [17] IRanges_2.18.3 S4Vectors_0.22.1
## [19] BiocParallel_1.18.1 Biobase_2.44.0
## [21] BiocGenerics_0.30.0 awst_0.0.4
## [23] readr_1.3.1
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 bit64_0.9-7 RColorBrewer_1.1-2
## [4] progress_1.2.2 httr_1.4.1 tools_3.6.3
## [7] R6_2.4.1 DBI_1.1.0 colorspace_1.4-1
## [10] gridExtra_2.3 tidyselect_1.0.0 prettyunits_1.1.1
## [13] bit_1.1-15.2 compiler_3.6.3 labeling_0.3
## [16] rtracklayer_1.44.4 scales_1.1.0 genefilter_1.66.0
## [19] askpass_1.1 DESeq_1.36.0 stringr_1.4.0
## [22] digest_0.6.25 rmarkdown_2.1 R.utils_2.9.2
## [25] jpeg_0.1-8.1 pkgconfig_2.0.3 htmltools_0.4.0
## [28] rlang_0.4.5 rstudioapi_0.11 RSQLite_2.2.0
## [31] farver_2.0.3 hwriter_1.3.2 jsonlite_1.6.1
## [34] dplyr_0.8.5 R.oo_1.23.0 RCurl_1.98-1.1
## [37] magrittr_1.5 GenomeInfoDbData_1.2.1 Matrix_1.2-18
## [40] Rcpp_1.0.3 munsell_0.5.0 viridis_0.5.1
## [43] reticulate_1.14 lifecycle_0.2.0 R.methodsS3_1.8.0
## [46] stringi_1.4.6 yaml_2.2.1 zlibbioc_1.30.0
## [49] plyr_1.8.6 grid_3.6.3 blob_1.2.1
## [52] crayon_1.3.4 lattice_0.20-41 splines_3.6.3
## [55] GenomicFeatures_1.36.4 annotate_1.62.0 hms_0.5.3
## [58] polylabelr_0.2.0 knitr_1.28 pillar_1.4.3
## [61] geneplotter_1.62.0 biomaRt_2.40.5 XML_3.99-0.3
## [64] glue_1.3.2 evaluate_0.14 latticeExtra_0.6-29
## [67] png_0.1-7 vctrs_0.2.4 polyclip_1.10-0
## [70] gtable_0.3.0 openssl_1.4.1 purrr_0.3.3
## [73] assertthat_0.2.1 ggplot2_3.3.0 xfun_0.12
## [76] aroma.light_3.14.0 xtable_1.8-4 RSpectra_0.16-0
## [79] viridisLite_0.3.0 survival_3.1-12 tibble_2.1.3
## [82] AnnotationDbi_1.46.1 memoise_1.1.0